Dynamic range of hypercubic stochastic excitable media 
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We study the response properties of d-dimensional hypercubic excitable networks to a stochas- 
tic stimulus. Each site, modelled either by a three-state stochastic susceptible-infected-recovered- 
susceptible system or by the probabilistic Greenberg-Hastings cellular automaton, is continuously 
and independently stimulated by an external Poisson rate h. The response function (mean density of 
active sites p versus h) is obtained via simulations (for d — 1, 2, 3, 4) and mean-field approximations 
at the single-site and pair levels (V d). In any dimension, the dynamic range and sensitivity of the 
response function are maximized precisely at the nonequilibrium phase transition to self-sustained 
activity, in agreement with a reasoning recently proposed. Moreover, the maximum dynamic range 
attained at a given dimension d is a decreasing function of d. 

PACS numbers: 87.19.L-, 87.10.-e, 87.18.Sn, 05.45.-a 



I. INTRODUCTION 

The building blocks of sensory organs are excitable 
neurons, upon which physical stimuli impinge continu- 
ously. Information about these stimuli is transformed 
into electrical activity of the neuronal membranes, usu- 
ally in the form of nonlinear excitations called spikes. 
The biophysics of ion channels involved in this process 
has been thoroughly investigated in the last five decades, 
after the success of the Hodgkin-Huxley theory [lj. De- 
spite immense progress in this regard, however, some fun- 
damental questions have remained unanswered. One of 
them has to do with two apparently conflicting experi- 
mental results. On the one hand, animals are subjected 
to stimulus intensities spanning many orders of mag- 
nitude, which their brains somehow manage to handle. 
This result is perhaps most easily revealed in psychophys- 
ical experiments. When humans are asked to assign an 
arbitrary (psychological) value to a given physical stimu- 
lus, this value is shown to be proportional to a power m of 
the stimulus intensity [1] (Stevens' law of psychophysics). 
The fact that the Stevens' exponent m is usually < 1 is 
consistent with the large dynamic range and sensitivity of 
psychophysical response functions (for instance, m ~ 0.6 
for the smell of heptane [5]). On the other hand, the 
response (mean firing rate) of isolated sensory neurons 
as a function of stimulus intensity has been shown to 
be an approximately linear saturating curve, at least for 
some sensory modalities. This implies that their dynamic 
range is usually small (for olfactory sensory neurons, they 
stay in the range of ~ 10 dB @,|j|). 

How is it then that large dynamic ranges are obtained 
from elements which individually have small dynamic 
ranges? What is the mechanism that generates Stevens' 
exponents m < 1? Two main mechanisms have been his- 
torically recognized as contributing to the phenomenon: 
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one of them invokes the intrinsic variation of thresholds 
in a population of sensory neurons The second one 
is adaptation, by which neurons adjust their operating 
ranges according to the statistics of the ambient stim- 
ulus (see, e.g., Refs. Jf| 0, H, 0, E3] for the case of the 
visual system). Both mechanisms certainly contribute 
to an enhancement of dynamic range. However, recent 
experimental data strongly suggest that additional mech- 
anisms based on a collective neuronal phenomenon could 
be at play: Deans et al. [ll| showed that knocking out 
gap junctions (electrical synapses among neurons) leads 
to a substantial change in the response function of retinal 
ganglion cells of mammals, with a decrease in dynamic 
range and sensitivity fllT ] and an increase in the response 
exponent 

In the last few years we have investigated this third 
possible mechanism, which addresses how neurons could 
cooperatively lead to an enhancement of dynamic range 
owing to the presence of lateral interactions (e.g. via 
chemical or electrical synapses). Connected, they form 
an extended system in which excitable waves are created 
upon incidence of incoming stimuli and annihilated upon 
collision (either with one another or with boundaries) due 
to the nonlinearity of their dynamics. The overall effect 
of this process is to collectively produce an enhancement 
of dynamic range and sensitivity, as compared to those 
of the elements alone [H [H M, M, QJ, Qj, [l| . We em- 
phasize that our proposal relies on very basic properties 
of excitable media, which opens the possibility that it 
could be applied not only to sensory systems, but wher- 
ever else enhanced sensitivity and dynamic range are re- 
quired. For example, extreme sensitivity is observed in 
rat motor cortex, where stimulation of a single pyramidal 
cell can evoke whisker movements [llj]. Also potentially 
related to what we propose is the experimental observa- 
tion that electrical synapses in the neocortex are present 
exclusively between inhibitory interneurons [2(| • Electri- 
cal coupling allows them to excite each other [20l ] and 
might have the functional role of augmenting their dy- 
namic range and sensitivity. 

The mechanism we discuss here is simple: let the ar- 
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rival of a suprathreshold stimulus reaching a sensory 
neuron be modelled by a Poisson process with rate h, 
which would be proportional to the stimulus intensity 
(say, the concentration of an odorant reaching the ol- 
factory epithelium). For very small h, stimulus events 
are rare and each of them would produce on average one 
excitation, if the excitable elements were disconnected. 
If they are connected, however, excitations can prop- 
agate stochastically to neighbors at some rate A. In 
this case, a single-stimulus event will generate an am- 
plifying excitable wave. If A is small, this wave will 
die out after some time, but the average network activ- 
ity p will nonetheless be larger than that of uncoupled 
neurons, leading to an enhanced sensitivity and larger 
dynamic range. This amplifying effect becomes more 
pronounced as A increases, so the dynamic range ini- 
tially increases with A. Increasing A further, however, 
one may reach a phase transition at some critical value 
A = A c , above which self-sustained activity becomes sta- 
ble [p(h = 0; A > A c ) > 0]. In this supercritical regime, 
the larger the coupling A, the more difficult it becomes 
for p to code for weak stimuli, which can hardly be dis- 
tinguished from the self-sustained background activity of 
the network. Therefore, for A > A c the dynamic range 
decreases with increasing A. Putting those two results 
together, one concludes that the dynamic range is maxi- 
mum at criticality (lrij . 

Clearly, the above reasoning applies to essentially any 
network topology. In its original version 16], it was for- 
mulated for an Erdos-Renyi random graph, where a sim- 
ple mean-field model perfectly captured the phenomenon. 
However, the enhancement of the dynamic range in that 
topology was about 50%, which is much less than what 
is observed experimentally (for instance, dynamic ranges 
in the olfactory glomerulus are at least twice as large 
as in olfactory sensory neurons [U H^]). This raises 
the question whether networks with different topologies 
could yield larger dynamic ranges and, if so, how these 
depend on network structural parameters. In Ref. [l8l |. 
for instance, the dynamic range of scale-free networks 
was shown to depend on the density of loops, but the 
evidence relied only on numerical simulations. In this 
contribution, we make use of simulations and analytical 
methods to deal with this question in hypercubic lattices, 
looking at the dimension d as a parameter. 

II. MODEL AND SIMULATION RESULTS 

We explore these ideas within a simple stochastic 
model of pulse-coupled excitable elements: 

P t (S x ) = -hP t (S :c ) - \J2Pt(S x ,Iy) + lPt(R x )(l) 

y 

Pt(4) = -P t (I x )+Xj2Pt(S x ,I y ) + hP t (S x ) (2) 
y 

P t {R x ) = -jPtiRx) + P t (h) , (3) 



where Pt(a x ) is the probability that site at location 
x G {1, . . . , N} is in state a at time t; Pt(a x , (3 V ) is the 
joint probability that sites at locations x and y are respec- 
tively in states a and (3 at time t; a, (3 G {5, /, R} denote 
a quiescent, excited or refractory state, respectively; y 
runs over the neighborhood of x; and 7 -1 is the charac- 
teristic refractory time, measured in units of the charac- 
teristic excitation time (defined as 1 , without loss of gen- 
erality [23]). We employ the notation of the stochastic 
susceptible-infected-recovered-susceptible (SIRS) model, 
to which this model is identical except for the external- 
stimulus field h, which is often missing in epidemiological 
modelling (where it amounts to a spontaneous infection 
rate [24]]). We can therefore extend a previous analysis of 
the stochastic SIRS model on a hypercubic lattice by Joo 
and Lebowitz [23j , from which our results can be derived 
if the external field h is added as in Eqs. dTJ- ((3]) . 

We are interested in the response function (or transfer 
function) of the excitable medium — i.e. the dependence 
of the stationary density of active sites p = lim^oo Pt(I) 
on the external stimulus intensity h (note that, in the 
context of neuroscience, the mean firing rate can be ob- 
tained by dividing p by the mean excitation time). Fig- 
ure [1] shows simulation results which confirm the general 
scenario described above. For 7=1 (which is kept fixed 
throughout this paper) and d = 2, a phase transition 
occurs at A = A c — 0.567(2). In the subcritical regime 
(A < A c ) p(h) is a linear-saturating curve whose slope in- 
creases monotonically with A. In the supercritical regime 
(A > A c ) one has po = lim?i— >o p{h) > 0, with po increas- 
ing monotonically with A [23j . 

To obtain the dynamic range of the response curve, we 
employ the definition usually adopted in the biological 
literature p. 1 . Let ft.0.1 and /iq.9 be the stimulus inten- 
sities such that their corresponding responses (po.i and 
P0.9) are, respectively, 10% and 90% above the base-line 
activity within the range [po , p m ax] '■ 

P-n = PO + V(pmax - Po) , (4) 

where p(h v ) = p v , p max = lim/^oc p(h) = 7/(7+1), and 
T) G {0.1, 0.9}. The dynamic range A is defined as 

asiobs '°0s)- < 5 » 

As depicted in Fig.QJa), A amounts to the range of stim- 
ulus intensities (measured in dB) that can be appropri- 
ately coded by the average activity in the network, dis- 
carding stimuli whose responses are too close either to 
baseline activity (p < po.i) or to saturation (p > po.g). 

According to this standard definition, the dynamic 
range of the response curves in Fig. [Ha) shows the pre- 
dicted behavior: for A < A c (A > A c ), A (A) is a mono- 
tonically increasing (decreasing) function. As shown in 
Fig. [21 the maximum dynamic range occurs precisely at 
criticality, where the response function is governed by 

r — 1 

the scaling relation p ~ h *■ . In all our simulations, the 
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FIG. 1: Response curves p(h) of the stochastic SIRS model 
on a two-dimensional lattice (results from simulations with 
N = 100 2 sites, periodic boundary conditions, averaged over 
a maximum time T max ~ 10 4 — 10 6 and five runs) in (a) 
linear-log scale and (b) log-log scale. From bottom to top, 
triangles denote A = 0, 0.2, . . . , 1. Filled circles denote A ~ 
A c . Relevant parameters for calculating the dynamic range 
A are exemplified in (a) for A = 1 (filled triangles), (b) Dot- 
dashed lines show literature values for the critical exponent 
6^ . Inset: response function near criticality for d — 1,3 and 4 
(system sizes and approximate critical coupling are N = 5000 



and A c ~ 7.73(8), N = 20 J and A c 
A c ~ 0.167(2), respectively). 



0.259(3), N = 10 4 and 



observed critical exponent S h is compatible with the 
literature values for the directed percolation (DP) uni- 



versality class: namely, S h = 0.111, 0.285, 0.45, and 
1/2 for d = 1, 2, 3, and d > 4, respectively [H, HE H3 
[see inset of Fig. HJb)] . This should be expected, since 
the model has local rules, a continuous transition to 
a unique absor bing state, and no further conservation 
laws [24 



, |29|, |30j . Because the exponent 5^ increases 
with increasing d and since it is the key element governing 
the dynamic range at criticality, we come to one of the 
main results of this paper: the maximum dynamic range 
attained at a given dimension d is a decreasing function 
of d. This result is summarized in Fig. [31 which exhibits 
the peaks of Fig. [5] versus the dimension of the lattice. 



We further note that similar results are obtained if 
one employs the Greenberg-Hastings cellular automaton 
(GHCA), where now the state transition S — > I is con- 
trolled by probabilities ph (external stimulus) and p\ 
(coupling), I — » R by p$ and R — > S by p 7 [311 ] . The 
only difference occurs for the particular case of excita- 
tions with a deterministic duration (ps = 1) m d = 1, in 
which case self-sustained activity cannot be stable and 
the maximum dynamic range occurs for p\ — 1 with 
an anomalous response exponent 1/2, as previously re- 
ported [IH (see stars in Fig. [3]). Note that ps = 1 seems 
biologically more realistic for neurons: while coupling 
may be stochastic, the duration of a spike is generally 
well described by a deterministic dynamics. For ps < 1 
and d = 1 a phase transition can occur just like in the 
SIRS model, and in this case A(A C ) is again a monotoni- 
cally decreasing function of d (as exemplified for ps = 0.5; 
see squares in Fig. |3|). 



III. 



MEAN-FIELD RESULTS 



A. Single-site approximation 

The remainder of this paper focuses on the possibili- 
ties of understanding these results with analytical means. 
Clearly, solving Eqs. HI)-© for A ^ is difficult, because 
the dynamics of single-site probabilities Pt(a x ) depend 
on two-site joint probabilities Pt(a x , y ), which in their 
turn depend on higher-order terms, and so forth. The 
simplest way to truncate this infinite hierarchy of equa- 
tions is the single-site (IS) mean field approximation, in 
which correlations are neglected [23j,[24(. The conditional 

15 

probabilities are approximated as Pt(a x \/3 y ) ~ Pt(oe x ), 
which implies that m-site joint probabilities get factor- 

ized as P(a x \\. . . , c&$) ~ 117=1 Pt{<&]), where a«> G 
{S, I, R}. In this approximation, assuming homogeneity 
and isotropy, Eqs. (P)-© reduce to a closed system and 
one obtains 



P(h) = S (^){ Pm \ x h + l-o) 




^OPmaxh 



{Pmaxh + 1 - a) 2 



(6) 



Note that information about the dimension d of the net- 
work is absorbed into an effective branching parameter 
a = Xz, where z = 2d is the number of neighbors each 
site has. When a — > 0, we obtain a linear saturating 
response p(h) = p m axh/ (p ma x + h), which is an exact 
result for uncoupled excitable elements. If h — > 0, Joo 
and Lebowitz's IS results are recovered: without an ex- 
ternal stimulus, the IS approximation predicts a phase 
transition at a = a c = 1 (A = A c = l/z), above which an 
active phase with po — p m axcr(cr — 1) > is stable (23[. 
The weak stimulus response is linear in the subcritical 

15 

regime, p(h; A < A c ) ~ h/{\ — zX). At criticality, how- 
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FIG. 2: Dynamic range vs. coupling: simulation results (symbols), IS (thin dashed) and 2S (thick dashed) mean field 
approximations. The peaks occur always at the phase transition. Half-widths of tuning of the critical regime [with heights 
measured relative to A(A = 0)] are AA ~ 1.4, 0.094, 0.059 and 0.042 for d = 1, 2, 3 and 4, respectively. 



IS 



1 ■ ' 1/2 

ever, it is governed by p(h; A = A c ) ~ (p m axh) , which 
leads to the mean-field exponent S^ 1 = 1/2. Applying 
the definition of Eqs. (0} and ([5]) to Eq. ©, one obtains 
the dynamic range in the IS approximation: 



A(a) 



is 



10 log 10 
101og 10 




(7) 



As depicted in Fig. [5J the peak at a = 1 is A(er c ) = 
301og 10 9 dB, which corresponds to an exact 50% en- 
hancement as compared to uncoupled elements: A(0) = 
201og 10 9 dB. Not surprisingly, this result holds for all d. 



B. Pair approximation 

To gain analytical insight into the dependence of A 
on d, we have solved Eqs. (P)-© in the so-called pair 
or two-site (2S) approximation, in which conditional 
probabilities are truncated beyond nearest neighbors: 

25 

P(a Xl \P X2 , Xx 3 ) ^ P(a Xl \(3 X2 ) (this approximation is 
valid for the hypercubic lattice because X3 is not a 
nearest neighbor of x±). This leads to three-site joint 

2S 

probabilities being approximated as P(a Xl , f3 X2 ,Xx 3 ) — 
P(a Xl , (3 X2 )P({3x2 ,Xx 3 )/ P(l3x 2 ) ■ To arrive at a closed set 



of equations, first one has to write down the dynam- 
ics also for two-site probabilities, which in this case are 
a direct extension of Joo and Lebowitz's equations for 
h 4- [2^. These equations ([ATT jl -lfAlB j) can be found 
in the Appendix. Applying the 2S approximation to all 
equations, one concludes, after some manipulation, that 
p satisfies a cubic equation [Eq. (|A.9[) ] . For h = one can 
obtain the 2S value of the coupling at which the phase 

transition occurs: A c = (7+I)/ (2d - 2 + (2d - 1)7) 

Figure [5] shows the dynamic ranges calculated from the 
numerically obtained response curves. The 2S approxi- 
mation shows a better agreement with simulations than 
does the IS approximation, and Fig. [5] clearly shows how 
the agreement improves with increasing d. In particular, 
note that the 2S approximation reproduces the inflection 
in the A(A) curve, which appears only for d — 1 [inset of 
Fig-I21 a )]- Like in the IS mean-field approximation and 
in the simulations, the subcritical response function for 
weak stimuli is linear: 



p(h;X< A c ) 



2S 



7 + 1 + A(7 + 2) 



[2d - 2 + (2d - 1)7] (A c - A) 



h. (8) 



More important for the topic of this article, the 2S ap- 
proximation manages to capture the dependence of A(A C ) 
on d, despite the fact that the weak-stimulus response at 
the critical coupling is still governed by a mean-field ex- 
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FIG. 3: Maximum dynamic range vs. dimension of the lat- 
tice: simulation results for the SIRS (open circles) and GHCA 
(stars denote ps = 1 and squares denote p$ = 0.5) models, IS 
(thin dashed) and 2S (thick dashed) approximations. GHCA 
simulations were performed with p 1 = 1/3 and the same num- 
ber of sites, Tmax and number of runs as described in the 
caption of Fig.[T]for the SIRS model. Inset: K(d) approaches 
the IS coefficient Pmax in the limit d — » oo (see eq. |9}. 
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FIG. 4: Active sites vs. time in a single run of a small d = 1 
SIRS system (one time step is counted every N = 200 up- 
dates; boundary conditions are open). An external stimulus 
with rate h = 10" 2 is applied for < t < 2000 (h = oth- 
erwise). Upper, middle and lower panels show the subcritical 
(A = 7.0), critical (A = 7.73) and supercritical (A = 8.5) 
regimes, respectively. 



ponent: p(\ c ;d) ~ K(d)h 1 ^ 2 . In this case, the decreasing 
function 



K(d) 



4d 2 7 ( 7 + l) 2 (2d- l)- 1 



7 3 (2d - 1) + 7 2 (6tf - 4) + 7 (6d - 5) + 2d - 1 

(9) 

correctly incorporates the influence of d on A(A C ) (see 
inset of Fig.©. 

The fact that neither the IS nor 2S approximations 
can reproduce the correct critical exponents for d < 4 



is well known. Inherent to those mean-field approxima- 
tions is the truncation of correlations, which is clearly 
inconsistent with the observed divergence of the correla- 
tion length £ and relaxation time r as criticality is ap- 
proached: £ ~ |A — A c | _!/ - L and r ~ |A — A r |~ t / H , where 



> and vu > are critical exponents 



In our 

simulations, correlation lengths and relaxation times at 
criticality are limited only by system size. At critical- 
ity also the survival probability V(t) decays as a power 
law (as opposed to an exponential fall in the subcriti- 
cal regime) [24j, which can lead to long-lived excitation 
waves, as illustrated in the single run of Fig. [U We note 
that correlations among cortical neurons several synapses 
distant from one another have been experimentally ob- 
served 32j and associated with the propagation of elec- 
trical waves [331 . 



IV. CONCLUDING REMARKS 

We have presented an analysis (with simulations and 
analytical results) of the response function of the stochas- 
tic SIRS model on hypercubic lattices. We confirmed 
that, as argued in Ref. [161 ]. the maximum dynamic range 
is obtained precisely at the nonequilibrium phase transi- 
tion where self-sustained activity becomes stable. More- 
over, since the response function at criticality is governed 
by the critical exponent S^ 1 , which for the DP universal- 
ity class increases with d, the maximum dynamic range 
obtained at a given dimension is a decreasing function 
of d. We therefore corroborate the claim that networks 
with spatial organization may have larger dynamic ranges 
than the random networks for which these ideas were first 
developed fill ]. 

This suggests the usefulness of low-dimensional arrays 
of excitable units for artificial sensor design, as well as 
raises speculations regarding the effective dimensional- 
ity of living neural networks. If one admits that large 
dynamic ranges could be favored by natural selection, 
organisms would tend to have their brains tuned at crit- 
icality, and in this case the mystery of how Stevens' ex- 
ponents < 1 arise would be solved: they would just be 
the critical exponent S^ 1 [HI]. 

In this context, our theoretical results join a recent 
flow of experimental evidence which is compatible with 
neurons collectively operating in a critical regime [Til [34l . 
IH, [H, H3, HI] . It is interesting to note that these exper- 
iments very often reveal exponents close to mean-field 
values: for instance, in cultures and acute slices of rat 
cortex, spontaneous activity occurs in avalanches whose 
size distribution decays as a power law with an exponent 
3/2 [34], HH, Hi| , which is the mean field result for branch- 
ing processes. Also, the response function of retinal gan- 
glion cells is well fitted by p ~ h 58 [llKIIj], yielding an 
exponent which is remarkably close to the Stevens' ex- 
ponent for light intensity and the DP S^ 1 exponent 
for d — 4. It would be interesting to investigate whether 
a small-world connectivity could conciliate these results, 
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on the one hand preserving the local order observed in 
real neural networks (see [39( for a recent review) while 
on the other hand allowing for mean-field exponents ow- 
ing to a small density of long-range connections. 
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Pt(S) = 1 -{ 1 + h)P t {S)- 1 P t {I) 

-z\P t (S,I), (A.4) 
P t (I) = hP t (S)-P t (I)+zXP t (S,I), (A.5) 
P t (S,I) = hP t (S)-(X + l + 2h)P t {S,I)-hP t (S,R) 

+(z-l)X^^[P t (S)-2P t (S,I) 

-Pt(S,R)}, (A.6) 
Pt (S, R) = 7 - jPt (S) - jPt (I) + P t (S, I) 
-(2<y +h)P t (S, R)-jP t (R, I) 
^ P t (S,I)P t (S,R) 

~ {z ~ 1)x — iw) — ' (AJ) 



APPENDIX: EQUATIONS FOR THE TWO-SITE 
APPROXIMATION 

As mentioned previously, the equations in Ref . (23j can 
be easily extended by including an external field h as 
follows: 

P t (S X ,Iy) = 1 P t {R X ,Iy)-{\+\ + h)P t {S X ,Iy) 

+ hP t (S X ,Sy)+ ^Pt(S X ,Sy,I W ) 

- *Pt{l W ,S X ,Iy) (A.l) 

Pt(S x ,Ry) = P t {S x ,I y ) + 1 P t {R x ,R y ) 

-{j + h)P t (S X ,Ry) 

XP t (I w ,S x ,Ry) (A.2) 

Pt{Rx,Iy) = -( 1 +l)P t (R X ,Iy)+P t (l X ,Iy) 
+ hP t (R x ,Sy) 

+ ^Pt(Rx,S y ,I w ), (A.3) 

where J\f x (y) is the neighborhood of y, excluding x. 
Note that we can omit the equation for P t (I x ,I y ) be- 
cause of the normalization condition Pt{A x , B y ) = 
Pt(By) [the same reasoning applies to P t (S x ,S y ) and 
Pt{R x , Ry)\- By applying the two-site approximation to 
Eqs. d)-© and Eqs. (|A.lj) - (|A.3j) under homogeneity and 
isotropy assumptions [23] , as well as the normalization 
condition P t (A x ) = 1, we obtain 



Pt{R,I) = P t (I)-P t (S,I) + hP t (S,R)-(2 + 7 )P t (R,I) 

The above equations form a closed system of ordi- 
nary differential equations. In its fixed point, p = 
lim^oc J* (J) is shown to satisfy the cubic equation 

A lP 3 + A 2 p 2 + A 3P + A A = , (A.9) 

where 

M = 7 2 { 7 3 [z 2 (z-l)A-z]+ 7 2 [z(2z 2 -2z-l)A-2z 

-1] + 7 [2z(z 2 - z - 1)A - 2z - 1] + z[(z 2 -z-l)X 
-1]} - /i{/i 2 ( 7 + I) 4 + htfz + 7 4 (zA + 4z + 3) + 

7 3 (4zA + 7z + 9) + 7 2 (6zA + 7z + 9) + 7 (4zA + 
4z + 3) + zX + z] + 7{7 4 (z 2 A + 2z) + 7 3 [zA(3z + 
2) + 6z + 3] + 2 7 2 [z(2z + 3)A + 4z + 3] + 3 7 [z(z + 

2)A + 2z+ 1] + z[(z + 2)A + 2)]}} , (A.10) 

A 2 = 7 jz 7 2 { 7 2 [-2z(z-l)A + z + l]+ 7 [-(3z 2 -4z 

-1)A + z + 3] - (2z 2 - 3z - 1)A + z + 1} + 

/i|3/i 2 ( 7 + l) 3 + /i[3 7 4 z + 7 3 (3zA + llz + 6) + 

7 2 (9zA + 16z + 12) + 7 (9zA + llz + 6) + 3zA + 
3z] + 7 {7 3 [3z 2 A + z(z + 4)] + 7 2 [4z(2z + 1)A + 
2z(z + 6) + 3] + 7 [z(9z + 8)A + 2z(z + 6) + 3] + 

z[4(z + l)A + z + 4]}}|, (A.ll) 
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A 3 = yj*V{7[(*-l)A-l] + («-2)A-l} 

1 A 4 = / l7 3 {/ l 2 ( 7 + l) + / l z[ 7 2 + 7 (A + 3) + A + l] + 

-/i{3^ 2 ( 7 + l) 2 + M3 7 3 z + 7 2 (3zA+10z + 3) + z 2 7 [ 7 (A + 1) + 2A + 1]}. (A.13) 

7 (6zA + lOz + 3) + 3zA + 3z] + z 7 { 7 2 [3zA + 

2(z + 1)] + 7 [(7z + 2)A + 3(z + 2)] + (5z + 2)A + 

^1 Cardan's formula yields the solution of eq. (|A.9j) . from 

2(z + 1)} > > , (A. 12) which the dynamic range can be numerically obtained. 
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